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Abstract. We present a new algorithm, called Multires- 
ohition Regularized Expectation Maximization (MREM) , 
for the reconstruction of 7-ray intensity maps from 
COMPTEL data. The algorithm is based on the itera- 
tive Richardson-Lucy scheme to which we added a wavelet 
thresholding step in order to eliminate image-noise in 
the reconstruction. The wavelet thresholding explicitly ac- 
counts for spatial correlations in the data, and adapts the 
angular resolution locally, depending on the significance 
of the signal in the data. 

We compare the performance of MREM to that of the 
maximum entropy and the Richardson-Lucy algorithms 
by means of Monte-Carlo simulations of COMPTEL 1.809 
MeV 7-ray line observations. The simulations demonstrate 
that the maximum entropy and Richardson-Lucy algo- 
rithms provide virtually identical reconstructions which 
are heavily disturbed by image noise. MREM largely sup- 
presses this noise in the reconstructions, showing only the 
significant structures that are present in the data. 

Application of MREM to COMPTEL 1.8 MeV 7-ray 
line data results in a 1.809 MeV sky map that is much 
smoother than the maximum entropy or Richardson-Lucy 
reconstructions presented previously. The essential fea- 
tures of this map are (1) an asymmetric galactic ridge 
emission reaching from I w 45° to I w 240°, (2) a bright 
localised emission feature in the Cygnus region around 
{l,b) w (80°, 0°), (3) two emission spots at I = 317° and 
I = 332° situated in the galactic plane, and (4) an ex- 
tended emission region around (l,b) « (160°, 0°). Compar- 
ison of the MREM map to the simulated reconstructions 
demonstrates that the 1.809 MeV emission is confined to 
the galactic plane. 
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1. Introduction 

Imaging the 7-ray sky at MeV energies by the COMPTEL 
telescope aboard the Compton Gamma-Ray Observatory 
(CGRO) presents a major methodological challenge. Reg- 
istered events are dominated by instrumental background, 
and additionally, source signals are widespread over the 
event parameter space. Consequently, image recovery re- 
lies on a complex deconvolution procedure and on the ac- 
curate modelling of the instrumental background compo- 
nent. A maximum entropy algorithm has been employed 
extensively for the reconstruction of intensity maps from 
COMPTEL data (Strong et al. 1992). Recent examples of 
maximum entropy all-sky maps can be found in Strong 
et al. (1999) and Bloemen et al. (1999a) for galactic con- 
tinuum emission or in Oberlack et al. (1996), Oberlack 
(1997) and Bloemen et al. (1999b) for 1.809 MeV 7-ray 
line radiation, attributed to the radioactive decay of 26 A1. 

Simulations revealed a tendency of clumpy reconstruc- 
tion of emission in our maximum entropy images, lead- 
ing to artificial 'hot spots' of 7-ray emission in the recon- 
structions of diffuse emission distributions (Knodlseder et 
al. 1996). From the images alone, these 'hot spots' are in- 
distinguishable from real point-like 7-ray sources, leading 
to considerable difficulties for the interpretation of the sky 
maps. Indeed, the assessment of the significance of indi- 
vidual 'hot spots' requires an substantial analysis effort, 
using simulations, Bootstrap analysis, and model fitting 
(e.g. Oberlack 1997). 

We understand the image lumpiness as the result of 
the weak constraints that are imposed on individual im- 
age pixels by our data. COMPTEL images are usually 
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reconstructed on a 1° x 1° pixel grid in order to exploit 
the telescope's angular location accuracy for point sources. 
The fine pixelisation implies, however, that for weak dif- 
fuse emission, 7-ray intensities in individual pixels are gen- 
erally not significant. Increasing the pixel size could in 
principle avoid this problem at the expense, however, of 
a reduced angular resolution. We want to notice that this 
is not a particular property of the maximum entropy al- 
gorithm, but of every method that is operating on a fixed 
grid of independent image pixels. Apparently, significance 
and angular resolution are intimately related quantities 
(this relation is more generally known as the bias-variance 
tradeoff) . 

Algorithms that rely on a pre-defined pixel grid re- 
quire an a priori choice of the angular resolution (by defin- 
ing the pixel size) without constraining the significance of 
the fluxes in individual image pixels. Alternatively, one 
may follow the opposite approach by choosing a priori the 
significance of image structures without constraining the 
angular resolution in the reconstruction. An implementa- 
tion of such an algorithm was discussed by Piha & Puet- 
ter (1993) who introduced generalised image cells (called 
'pixons') to correlate adjacent image pixels according to 
the signal strength. However, the application of Pixon 
based image reconstruction to COMPTEL 1.8 MeV did 
not provide satisfactory results (Knodlseder et al. 1996). 

In this paper we present a new algorithm, called 
Multiresolution Regularized Expectation Maximization 
(MREM), which we developed in particular for the re- 
construction of diffuse 7-ray emission. We combined an 
expectation maximization (EM) algorithm with a multi- 
resolution analysis based on wavelets, which explicitly ac- 
counts for spatial correlations in the reconstructed image. 
This leads to a convergent algorithm which automatically 
stops when the significant structure has been extracted 
from the data (by significant structure we mean struc- 
ture that will not change much under perturbation of the 
data). The method requires an a priori choice of the sig- 
nificance level of emission structures while adapting the 
angular resolution according to the signal. 

In the following we will present the MREM algorithm 
(^) and illustrate its performance by means of simulations 
of COMPTEL observations (§|). The MREM algorithm is 
then applied for the reconstruction of an 1.809 MeV all- 
sky map based on COMPTEL data obtained between May 
1991 to June 1996 (§[|). This sky map will be compared 
to 1.8 MeV all-sky maps presented previously which have 
been derived by the maximum entropy method (Oberlack 
et al. 1996) or the Richardson-Lucy algorithm (Knodlseder 
et al. 1996). A more theoretical description of the MREM 
algorithm will be given in a separate paper (Dixon et al., 
in preparation). 



2. The MREM algorithm 

MREM is based on the Richardson-Lucy (RL) algorithm 
which has been proposed by Richardson (1972) and Lucy 
(1974) for the restoration of degraded images. Given an 
initial estimate f® for the image, RL iteratively improves 
this estimate using 



fk+l _ rk 



EN m td 
i=i e k r H: 



(1) 



(k denotes the iteration). By is the instrumental response 
matrix which links the data space (indexed by i) to the 
image space (indexed by j). For a given image f k and a 
given background model hi, the expected number of counts 
in a data space cell is given by e\ = Ylj=i Rijfj + The 
number of events observed in data space cell i is given by 
rij. It is easily seen that Eq. (Q) may be also written in the 
additive form 
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where 
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is the additive RL correction. 

Shepp & Vardi (1982) demonstrated that the 
Richardson-Lucy scheme is a special case of the expecta- 
tion maximization (EM) algorithm (Dempster et al. 1977), 
and consequently it converges to the positively constrained 
maximum likelihood solution for Poisson data. Due to the 
slow convergence of the algorithm, several modifications 
have been proposed to accelerate convergence (Fessler & 
Hero 1994). For COMPTEL data we found that the ML- 
LINB-1 algorithm of Kaufman (1987) gives reasonable ac- 
celeration without degrading the reconstruction proper- 
ties. For ML-LINB-1, Eq. (§) is replaced by 

f- +1 = f- + A fe 5/j=, (4) 

where X k is determined for each iteration using a line- 
search in order to maximise the likelihood for fj +1 subject 
to the constraint X k 8f k > —fj (this constraint ensures the 
positivity of the intensities) . 

It is obvious from Eqs. (0) to (0) that RL operates on 
a pre-defined pixel grid without any direct correlation be- 
tween individual pixels. In particular, apart from the con- 
volution with the transpose of the response matrix By, 
there is nothing which prevents RL from matching the 
estimates e\ to the measurement rij, and noise can eas- 
ily propagate into the reconstruction where it is generally 
amplified. 

For these reasons we added a multiresolution analy- 
sis to the iterative procedure which aims to correlate the 
image pixels and to extract only the significant structure 
from the data. 
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Each iteration of our MREM algorithm is composed 
of four steps: First, we evaluate the normalised correction 
map 



5K 



Ef =1 (f - i) R, 



(5) 



for which Var(5h k ) = 1. Second, 5h k is transformed into 
the wavelet domain where it is represented by a set of 
wavelet coefficients w l m , I representing the scales, and m 
denoting the wavelet coefficients at this scale. At scales 
I > 1 (I = 1 represents a DC offset) the coefficients falling 
below a given threshold t 1 are zeroed by applying the 
operator 



r)(w,T l ) 







M < r 
\w\ > r 



(6) 



This method is generally referred to as wavelet threshold- 
ing and has been proven successful for the removing of 
noise from a dataset without smoothing out sharp struc- 
tures (Donoho 1993; Graps 1995). Backtransformation of 
the nonzero coefficients from the wavelet domain into the 
image domain provides then a de-noised correction map 
Sh k 
by 



; . In compact matrix notation, the second step is given 



Sh k = W T nW5h h 



(7) 



where W is the discrete wavelet transform (throughout 
this paper we use the translation invariant 'cycle spinning' 
transformation of Coifman & Donoho (1995) and employ 
Coiflet wavelets with 4 parameters). Third, we calculate 



6 ft = ftSh) 



V 



(8) 



which, in absence of any wavelet thresholding, is equiva- 
lent to the original RL correction map Eq. (j|) . In the last 
step, the previous estimate /j is updated using Eq. (|J). 
Due to the wavelet thresholding, positivity of the pixel in- 
tensities is not implicitly assured, and we explicitly require 
jk+i > where f e is a negligible intensity level. 

For efficient de-noising, the scale-dependent thresholds 
t have to be related to the expected statistical noise a 
in the wavelet domain at each scale. We estimate a 1 by 
simulations where we replace n; in Eq. (^]) by a Poisson 
derivate of e\ and transformation of the resulting 'mock 
correction map' into the wavelet domain. In this approach 
it is important that the statistical noise in the MREM 
correction map 5h k is independent of the pixel location j. 
For this reason we normalised 8h k so that Var(Shj ) = 1. 
We define then t 1 — so 1 , where s specifies the significance 
level below which structures should be suppressed in the 
reconstruction. In the examples presented below we will 



vary s between 2.5 and 3.5 in order to demonstrate the 
impact of the choice of s on the reconstructed images. 

The final critical aspect of MREM is that we do not 
calculate corrections for all wavelet scales simultaneously. 
We begin by allowing in only corrections corresponding 
to the largest wavelet scale, with the other scales simply 
being zeroed out; thus our estimate in the initial iterations 
corresponds only to the average large scale structure in the 
map. Once this converges, we then admit corrections from 
the next smallest wavelet scale, allow it to converge, and 
so forth. The resulting reconstruction is the final product 
of the algorithm, which generally requires between 20 and 
40 iterations for our data. 

This procedure of progressively admitting smaller 
wavelets is crucial to the performance of MREM, and 
has some rather interesting ramifications which we shall 
discuss in detail in a subsequent paper (Dixon et al., in 
preparation). We briefly note here that its main purpose 
is to aid in the discrimination of noise-induced corrections 
from that structure which is "stable" in the sense of be- 
ing reproducible from different datasets. Examination of 
Eq. (|^) indicates that if the estimates e\ are very differ- 
ent from the data n i; the corresponding corrections will 
also be large. Ideally, we would like this to occur only if 
the correction corresponds to some statistically interesting 
structure, but for arbitrary e\ (e.g., the initial flat guess) 
this won't be the case. By first converging to a coarse 
approximation, we get an estimate that is "close" to the 
next coarsest approximation, which generally forces the 
noise-induced corrections at that next smallest scale to be 
small compared to those which we deem "interesting" . It 
is further interesting to note that the unregularised RL 
iteration tends to pick out the larger scale average struc- 
ture first, only adding details in later iterations, and in 
this sense the progressive scale procedure dovetails nicely 
with the known characteristics of the iteration. 

The de-noising using the wavelet transform has the de- 
sired property of introducing pixel-to-pixel correlations in 
the image where the correlation length depends on the 
amount of structure in the data. Regions of the sky with 
uniform emission will be represented by few large-scale 
wavelet coefficients, while point sources are represented by 
few small-scale coefficients. An important feature of our 
algorithm is that it is convergent. If all significant struc- 
ture has been extracted from the data, where 'significant' 
is defined by the choice of s, the thresholding operator 
will zero all wavelet coefficients and consequently the cor- 
rection map will be structureless. At this point, further 
iterations won't alter the reconstructed image anymore, 
hence we stop the iterations. 

3. Simulations 

To illustrate the performance of MREM with respect to 
the maximum entropy (ME) and the Richardson-Lucy 
(RL) algorithms, we apply them to simulated COMP- 
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TEL observations of 1.809 MeV 7-ray line emission. The 
mock data that are used in the simulations are based on 
a two-component data space model, composed of the in- 
strumental background and adopted models for the 7-ray 
line distribution for two typical cases: a smooth large-scale 
emission model, and a rather structured model with emis- 
sion on many spatial scales. The instrumental response 
and background were calculated as expected for the com- 
bination of observation periods 0.1 — 522.5, correspond- 
ing to data taken between May 1991 to June 1996. From 
both components of the data space model mock datasets 
were created independently by means of a random num- 
ber generator assuming Poisson noise. Both components 
were then added and images have been reconstructed from 
the combined mock dataset. For the reconstructions it has 
been assumed that the instrumental background is known 
precisely, hence the resulting images are not subject to 
possible systematic uncertainties of the employed back- 
ground model. They are sensitive, however, to statistical 
uncertainties which are due to the particular data 'real- 
isation' as obtained by the random sampling procedure. 
To illustrate this sensitivity, the same mock dataset has 
been used for the instrumental background component in 
all simulations. 

The following 1.809 MeV model intensity distributions 
have been chosen. First, we use an exponential disk model, 
i.e. the intensity distribution that is expected if the galac- 
tic 26 A1 mass density would follow a double exponential 
law with scale radius of Rq =4.5 kpc and scale height of 
zq = 90 pc. Model fitting has confirmed that these param- 
eters provide a reasonable first-order description of 1.809 
MeV emission (Knodlseder 1997). The total galactic 26 A1 
mass has been normalised to 3Mq , a value slightly in ex- 
cess from recent findings (e.g. Diehl et al. 1997). Second, 
the EGRET > 100 MeV all-sky map was taken as tem- 
plate for the 1.809 MeV intensity distribution. The 1.809 
MeV intensity level of the map was adjusted to a plausible 
level by fitting the map to COMPTEL 1.8 MeV 7-ray line 
data. The first case testifies the response of the image re- 
construction algorithms to a smooth intensity distribution 
while the second case represents probably a more realistic 
situation with structure on all spatial scales, from point- 
like to diffuse galactic plane emission. 

3.1. Exponential disk model 

The results of the exponential disk simulation are com- 
piled in Fig. [l]. For comparison, the intensity distribution 
of the exponential disk model is also shown. Since in our 
implementation ME and RL provide no criteria where to 
stop the iterative procedure, we used the correlation co- 
efficient between the reconstruction and the model inten- 
sity map to determine the iteration which provides the 
smallest discrepancy to the model. This is the case after 
iteration 8 for ME and the accelerated RL algorithm. 



Both the ME and the RL reconstructions clearly 
pick up the emission ridge along the galactic plane with 
the highest intensities found towards the central radian 
(—30° < I < 30°). The most striking difference between 
the model and the ME and RL reconstructions, however, 
is the lumpiness of the recovered sky maps. Although the 
emission follows in average the model intensity profile, it 
exhibits strong oscillations around this average, leading 
to 'hot spots' and emission gaps along the galactic plane. 
Indeed, these oscillations are already present in the first 
iterations of the reconstruction process and become more 
and more amplified with proceeding iterations. If the it- 
erations are pursued beyond those shown in Fig. [l], the 
oscillations will break up, and the image will be composed 
of nearly isolated point sources (Knodlseder et al. 1996). 

It is also interesting to recognise that the ME recon- 
struction is virtually identical to the RL reconstruction. 
The difference between both algorithms is that ME im- 
poses an additional constraint on the reconstructed im- 
age in that it 'pushes' the image towards a 'flat' sky map 
- especially if the data are not very constraining. This 
results in systematically lower fluxes for the ME recon- 
structions with respect to RL, which can be seen from the 
intensity profiles in Fig. [l]. If the ME iterations are pro- 
ceeded further, and hence the entropy criterion is gradu- 
ally weakened, the flux discrepancy between ME and RL 
disappears. For this reason we always use 'high' (e.g. 20- 
30) iterations when we determine fluxes from our ME sky 
maps. 

In contrast to ME and RL, the MREM reconstruc- 
tions provide rather smooth emission distributions. While 
some lumpiness remains for s = 2.5, the images obtained 
with s = 3.0 and 3.5 show no 'hot spots' or emission 
gaps. In particular, the longitude profile is reasonably 
well reproduced and obeys only small deviations from the 
model distribution. The most striking difference between 
the MREM sky maps and the model is the larger latitude 
extent of the reconstructions. This, however, is not surpris- 
ing since the width of the exponential disk model of 2.7° 
(FWHM) is considerably smaller than the instrument's 
angular resolution of 4° (FWHM) at 1.8 MeV (Schonfelder 
et al. 1993). Together with the weakness of the signal, this 
limits the achievable resolution in the reconstructions. In- 
deed, the width of the latitude profile depends on the se- 
lected significance level s, rising from 5.3° for s = 2.5 
to 9.6° for s = 3.5. Obviously, the significance of the re- 
covered emission features and the angular resolution are 
intimately related quantities. Note that the width of the 
latitude profiles obtained by ME and RL is 5.7° (FWHM), 
which is also considerably wider than that of the model. 

To judge the quality of the reconstructed images we 
determine the 1.8 MeV 7-ray line residuals by means of 
a maximum likelihood ratio test (de Boer et al. 1992). 
For this purpose the sky maps of Fig. [I] are convolved into 
the COMPTEL data space and added to the instrumental 
background model. Residual emission is then searched by 



J. Knodlseder et al.: Image Reconstruction of COMPTEL 1.8 MeV 26 A1 Line Data 5 




180 150 120 90 60 30 -30 -60 -90 -120 -150 -180 

Galactic Longitude (deg) 




Galactic Latitude (deg) 

Fig. 1. Reconstruction of an exponential disk intensity distribution (Exp. -disk) using the maximum entropy method 
(ME), the Richardson-Lucy algorithm (RL), and the Multiresolution Regularized Expectation Maximization algorithm 
(MREM). MREM reconstructions for the significance levels s = 2.5, 3.0, and 3.5 are shown. The bottom plots compare 
the longitude profile and latitude profile of the reconstructions to that of the exponential disk model; the model is shown 
as solid, the ME reconstruction as dotted, the RL reconstruction as dashed, and the MREM (s = 3.0) reconstruction 
as dashed-dotted line. 
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Fig. 2. Residual maximum likelihood ratio maps of the exponential disk simulations. Contour levels: \/ — 2 In A = 
2,3,.. .. From top to bottom the panels show (a) the residuals of the background sample only (i.e. the noise), (b) the 
residuals of iteration 8 of the ME reconstruction, (c) the residuals of iteration 8 of the RL reconstruction, and (d) the 
residuals of the MREM (s = 3.0) reconstruction. 



fitting point source models on top of the combined data 
space model for a grid of source positions. The results of 
this point source search are shown in Fig. |^. The quantity 
plotted is —2 In A, where A is the maximum likelihood ra- 
tio L(M)/L(S + M), M represents the (two-component) 
data space model, and 5* the source model which is moved 
over the sky area searched for residual emission. In such 
a search, —2 In A obeys a xi distribution; in studies of a 
given source, xi applies. In the latter case, the point source 
significance (in Gaussian a) is given by y— 2 In A. 



The top panel in Fig. Q shows the residuals of the in- 
strumental background sample only, hence reflects the sta- 
tistical noise in the mock datasets (due to the dominance 
of the instrumental background component, the statistical 
noise is dominated by the background fluctuations). In the 
ideal case, the residuals of the reconstructed images should 
be almost identical to those of the background sample. In- 
deed, the residuals found on top of the MREM (s = 3.0) 
reconstruction (panel d) are very similar to those expected 
for an ideal reconstruction (panel a) . The features are ba- 
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sically identical; only small deviations are found in their 
amplitude, e.g. at I i=a 110° where MREM slightly over- 
estimates the emission. The ME reconstruction (panel b) 
forces image flatness, hence the residuals reflect a promi- 
nent flux suppression of the entire plane emission. We 
therefore cannot easily detect to which extent noise has 
been included in the ME reconstruction. For RL (panel 
c), no residual 1.8 MeV emission is seen that correlates 
with the galactic plane. In contrary, the likelihood ra- 
tios are even too small for the RL reconstruction with 
respect to the noise simulation, as expected if the data 
were overfit by the intensity map. Additionally, promi- 
nent background features, like those at (l,b) = (113°, 2°) 
or at (l,b) = (81°,— 16°), are drastically reduced in both 
the ME and RL residual maps, yet are perceptible in the 
reconstructed intensity maps (cf. Fig. [l). This illustrates 
that statistical noise in the data is at least partially fit by 
the ME and RL sky maps. It follows that the lumpiness 
of the ME and RL reconstructions are due to overfitting 
of the data. 



3.2. EGRET > 100 MeV sky map 

Figure ^ presents the results of the simulations based on 
the EGRET > 100 MeV all-sky map. This map shows a 
ridge of diffuse emission along the galactic plane with a 
notably intensity enhancement towards the inner Galaxy, 
some prominent galactic point sources, some localised 
emission regions, and some extragalactic point sources. 
The intensity level of the EGRET > 100 MeV all-sky map 
has been adjusted to a plausible 1.8 MeV intensity level 
by fitting the map to COMPTEL 1.8 MeV 7-ray line data, 
resulting in a 1.809 MeV line flux from the inner radian 



< 20°) of - 3 x 10~ 4 ph cm' 



L rad 



This ad- 



justment pushed most of the point sources in the EGRET 
map below the sensitivity limit of COMPTEL at 1.8 MeV, 
leaving Vela with 3 x 10~ 5 ph cm _2 s _1 and Geminga with 
1 x 10~ 5 ph cm _2 s _1 as the most prominent objects. 

Indeed, the only point source recovered in the ME and 
RL reconstructions is Vela (I — 264°, b = —3°), while 
only a small hint of 7-ray emission is seen at the posi- 
tion of Geminga (I = 195°, b = —4°). Similar to the ex- 
ponential disk simulation, 'hot spots' and emission gaps 
appear along the galactic plane which only occasionally 
coincide with localised features in the model map. Such 
coincidences are found e.g. at I « 80° (Cygnus) or at 
I « —45°. However, localised emission features are also 
found in the exponential disk simulations at these posi- 
tions (cf. Fig. |l|) where no such features are present in 
the model. Since only the mock dataset of the (dominant) 
instrumental background component is common to both 
simulations, it is very suggestive that the observed fea- 
tures are at least partially due to positive statistical fluc- 
tuations of the background data. Other localised emission 
features in the EGRET map, like the spots at I ~ 20°, 
I —18°, or I —75° (Carina), coincide with negative 



statistical fluctuations of the background component and 
annihilate; consequently no feature is seen in the recon- 
structions at these positions. Additionally, artificial 'hot 
spots' appear in the ME and RL maps where no such fea- 
tures are present in the model. Examples are the strong 
feature at I « 115° and the spur towards negative latitudes 
at I 1=3 —15°. Again, these features can also be percepted 
in the exponential disk simulations, confirming that they 
arise from the statistical noise of the background sample. 

In contrast to ME and RL, MREM again provides 
much smoother reconstructions of the data, avoiding most 
of the artifacts. In the s = 2.5 run, only the most promi- 
nent artifacts are visible (e.g. the spot at I « ff5°), but 
many of the real localised features are recovered (I « 80° , 
I pa 135°, I w —45°, and Vela). Increasing the requirement 
for the significance of the emission structures to s = 3.0 
removes the remaining artifacts, but eliminates also most 
of the localised emission features. Nevertheless, weak hints 
for Vela and the I « 135° source are still present in the 
sky map. These hints disappear when s is increased to 3.5. 
Again, the latitude extent of the reconstructions is slightly 
higher than that of the models due to the combined result 
of the instrument's angular resolution of only 4° (FWHM) 
together with a low signal to noise ratio. Yet, the extended 
diffuse emission above the galactic centre is still recovered 
in the maps. 

The residual analysis of the MREM s = 3.0 recon- 
struction reveals only weak emission at the position of the 
localised features, indicating that they are not very signif- 
icant (cf. Fig. ||). The most prominent residuals are found 
at the position of Vela and at (l,b) = (81°,— 16°), with 
likelihood ratios of —2 In A = 16.6 and 14.3, respectively. 
While the first residual corresponds to a real source in 
the EGRET map, the second one is a clear background 
fluctuation. If the existence of the Vela source would not 
be known a priori, the likelihood ratio of 16.6 converts to 
a detection significance of 3.3er (3 d.o.f.). Taking into ac- 
count the number of trials made in the point source search, 
this value can not be interpreted as a significant detection. 
However, if the Vela source is considered as known object, 
the likelihood ratio converts to a 4.1cr detection signifi- 
cance (1 d.o.f.). The major objective of an 1.809 MeV all- 
sky map, however, is the discovery of unknown objects, 
hence it is desirable that the Vela source is not recovered 
in the reconstruction. Otherwise, as demonstrated by the 
ME and RL or the MREM (s = 2.5) reconstructions, arti- 
facts will also enter the reconstruction, making the inter- 
pretation of the sky map difficult. 

To illustrate that MREM indeed recovers point sources 
if they are significant, we performed an additional simu- 
lation where we increased the intensity of the EGRET 
> 100 MeV template by a factor of 5 with respect to the 
1.809 MeV intensity. This corresponds to an increase of 
a factor of 5 in the signal-to-noise ratio, which is equiva- 
lent to a sensitivity enhancement of the same magnitude. 
The resulting MREM reconstruction is shown in Fig. || 
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Fig. 3. Reconstruction of an EGRET > 100 MeV template which was adjusted to the expected 1.8 MeV intensity 
by means of a model fit. In the profiles, the model is shown as solid, the ME reconstruction as dotted, the RL 
reconstruction as dashed, and the MREM (s = 3.0) reconstruction as dashed-dotted line. 
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Fig. 4. Residual maximum likelihood ratio maps of the EGRET > 100 MeV all-sky map simulations. Contour levels 
like in Fig. ||[ From top to bottom the panels show (a) the residuals of the background sample only (i.e. the noise), 
(b) the residuals of iteration 6 of the ME reconstruction, (c) the residuals of iteration 6 of the RL reconstruction, and 
(d) the residuals of the MREM (s = 3.0) reconstruction. 



for s — 3.0. As expected, much more structure along the 
galactic plane is now recovered. Prominent point sources, 
such as Vela, Geminga, or the Crab, and localised emission 
features, e.g. in Cygnus and Carina, are now clearly visi- 
ble. This demonstrates that the absence of these features 
in the MREM map derived for the 1.809 MeV intensity 
level (Fig. |) relates to their significance. 



4. The COMPTEL 1.8 MeV sky 

The MREM algorithm is now applied to real COMPTEL 
1.8 MeV data, taken during observation periods 0.1 - 522.5 
(May 1991 to June 1996). The instrumental background 
component was estimated using contemporaneous data 
at adjacent energies, following the procedure described 
in Knodlseder et al. (1996). The time- variability of the 
instrumental background was taken into account by de- 
termination of the background model on single observa- 
tion basis, and by its proper relative normalisation using 
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Fig. 5. MREM (s = 3.0) reconstruction of an EGRET 
> 100 MeV template with an intensity that is 5 times 
larger than the 1.809 MeV intensity. 

the activation history of major background components 
(Oberlack 1997). Since the absolute normalisation as well 
as the if distribution of the background model are only 
weakly constrained, we added an additional step where we 
determine both by all-sky model fitting (maximum likeli- 
hood optimisation). For this purpose we fitted the instru- 
mental background model together with a template for 
the 1.809 MeV intensity distribution to the COMPTEL 
1.8 MeV data, where we determined independent scaling 
factors for all (p layers of the background model as well as 
a global scaling factor for the 1.8 MeV intensity template. 
This procedure provides an estimate of the total 1.809 
MeV sky flux and an improved estimate of the instrumen- 
tal background component, which is then used for image 
reconstruction. For the 1.809 MeV template we used the 
53 GHz free-free emission map derived from COBE/DMR 
data (Bennett et al. 1992) which was found to provide 
the best description of the COMPTEL 1.8 MeV data in a 
recent study using a wide variety of models (Knodlseder 
et al. 1999). An alternative method for deriving an in- 
strumental background model for the analysis of 1.8 MeV 
7-ray line data is described in Bloemen et al. (1999b). 

Figure | shows the COMPTEL 1.809 MeV 7-ray line 
all-sky maps that are obtained by the different reconstruc- 
tion algorithms. The maximum entropy and Richardson- 
Lucy maps are similar to those presented in previous work 
(Diehl et al. 1995; Oberlack et al. 1996; Knodlseder et 
al. 1996; Oberlack 1997; Bloemen et al. 1999b) with minor 
differences being due to differences in the analysed data 
volume or the employed background modelling procedure. 
The most distinct feature in these maps is emission along 
the ridge of the galactic plane. Again we see the lumpiness 
that our above simulations also show for these methods, 
indicating overfitting of the data. According to the dis- 
cussion above we cannot decide from the sky maps alone 
which of the lumps may correspond to real emission and 
which are artifacts due to the background noise. 

MREM avoids this confusion, suppressing efficiently 
the noise components of the image. The reconstructed in- 
tensity profiles are characterised by a notable asymme- 



try with respect to the galactic centre and some localised 
emission features. The most prominent of these features 
is located in the Cygnus region around (l 7 b) » (80°, 0°) 
where a bright extended emission spot is clearly separated 
from the inner galactic ridge emission by a bridge of rela- 
tively low 1.8 MeV intensity. The same feature is also seen 
in the ME and RL maps where it obeys a much more com- 
plex structure. The MREM reconstructions suggest that 
most of this structure is not individually significant, and 
could as well be more diffuse or located differently. We 
therefore safely may extract the fact of significant Cygnus 
region emission, separated from the inner galactic ridge. 

This inner ridge appears smooth on the MREM im- 
age, yet also here reveals a pronounced asymmetry with 
respect to the galactic centre. While at positive longitudes 
the intensity drops steeply from I rj 30° to I w 50°, the 
1.8 MeV emission extends continuously to I 1=3 240° at 
negative longitudes. Along the ridge the MREM recon- 
structions reveal only little structure. From the many 'hot 
spots' seen in the ME and RL reconstructions, only the 
most prominent ones are still perceptible in the MREM 
image obtained for s — 2.5. Increasing the significance 
level to s — 3.0 removes most of them, keeping only two 
emission spots at I — 317° and 332° which are separated 
by a weak emission gap at I — 324°. This is the most per- 
sistent structure along the inner galactic plane ridge which 
is still clearly visible for a significance level of s = 3.5. It 
is also very pronounced in the ME and RL maps. Ad- 
ditional hints for weak excess emission are found in the 
longitude profile at I = 21°, 30°, 44°, 286°, and 345°, but 
they disappear if s is increased to 3.5. Obviously, the sig- 
nificance of these excesses is close to the sensitivity limit 
of COMPTEL, and the assessment of their reality needs 
more dedicated studies. 

The distinct emission gap which separates two lo- 
calised emission regions at I — 266° (Vela) and I — 286° 
(Carina) in the ME and RL maps is not seen in the MREM 
reconstructions. Yet, a weak intensity dip is found in the 
s = 2.5 and 3.0 MREM maps at this location, indicating 
that some structure may indeed be present. The promi- 
nent 'hot spot' towards the galactic centre, clearly visible 
in the ME and RL maps at I — 4°, is only present in the 
MREM map for s — 2.5, but disappears for higher signifi- 
cance levels. Apparently, the data are also consistent with 
a smooth emission profile in this region. 

Comparison of the 1.8 MeV maps with the EGRET 
map simulations indicates that the 1.809 MeV emission is 
very confined to the galactic plane. In particular, there 
is no hint for an extended emission component similar to 
that seen in the EGRET map above the inner Galaxy. In- 
deed, model fitting using exponential disk models revealed 
a small scale height of zq = 90 pc for the galactic 26 A1 dis- 
tribution (Knodlseder 1997). Comparison of the 1.8 MeV 
maps with the exponential disk simulations (for which a 
scale height of zq = 90 pc was assumed) confirms this re- 
sult. In particular, the width of the MREM (s = 3.0) 1.8 
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Fig. 6. COMPTEL 1.8 MeV all-sky maps derived using the maximum entropy method (ME), the Richardson-Lucy 
algorithm (RL), and the Multiresolution Expectation Maximization algorithm (MREM). In the profiles, the MREM 
(s = 3.0) reconstruction is shown as solid, the ME reconstruction as dotted, and the RL reconstruction as dashed line. 
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Fig. 7. Residual maximum likelihood ratio maps of the COMPTEL 1.8 MeV all-sky maps. Contour levels like in 
Fig. From top to bottom the panels show (a) the residuals of iteration 6 of the ME reconstruction, (b) the residuals 
of iteration 6 of the RL reconstruction, and (c) the residuals of the MREM (s = 3.0) reconstruction. 



MeV latitude profile (7.5° FWHM) is even smaller than 
that obtained for the exponential disk simulation (8.8° 
FWHM), indicating that the scale height of the 26 A1 dis- 
tribution may be even below 90 pc. 

Near the anticentre, all maps of Fig. ^ show indications 
for extended 1.8 MeV 7-ray line emission. In the ME and 
RL reconstructions, weak emission spots are spread over 
a region extending from 125° — 170° in galactic longitude 
and from -20° - 30° in galactic latitude. The MREM 
algorithm combines these spots to a more concentrated 
emission structure, roughly located at I ~ 160° with an 
angular extent of ~ 20°. This again illustrates that the 
spots in the ME and RL images are not significant for 
themselves, but when combined they provide a significant 
1.809 MeV emission feature. 

Residual maximum likelihood ratio maps of the 
COMPTEL 1.8 MeV all-sky maps are compiled in Fig. ^. 
The ME reconstruction shows significant residual emission 
along the galactic plane which is strongly correlated to 
the reconstructed sky intensity distribution. The intensity 
profiles in Fig. ^ illustrate that iteration 6 of the ME recon- 



struction considerable underestimates 1.8 MeV intensities 
with respect to RL and MREM. For higher ME iterations, 
this underestimation disappears as the maximum entropy 
reconstruction approaches the maximum likelihood solu- 
tion. Yet, the diffuse intensity distribution breaks up into 
nearly isolated point sources for late iterations due to over- 
fit of the data (Knodlseder et al. 1996). Therefore we typ- 
ically present COMPTEL ME images and longitude pro- 
files from 'early' iterations in order not to emphasise ar- 
tificial structures, while 'late' iterations are used to de- 
rive 1.809 MeV fluxes and latitude profiles in order to 
recover the correct flux values (Diehl et al. 1995; Ober- 
lack et al. 1996). Alternatively, intensity distributions for 
'late' iterations have been smoothed to the instrumental 
resolution for image presentation to reduce the artificial 
lumpiness (Oberlack 1997; Strong et al. 1999). 

Also the RL reconstruction shows residuals that are 
correlated with the galactic plane, although they are much 
smaller than for ME. Yet there are regions where almost 
no residuals are found, in particular at negative galactic 
longitudes (I < 0°) above and below the galactic plane. 
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Comparison with the simulations suggests that the lack of 
residuals is again due to overfit of the data. In contrast, the 
MREM (s = 3.0) reconstruction provides residuals that 
appear uncorrelated with the galactic plane. This clearly 
illustrates that the MREM map presents a statistical sat- 
isfactory description of COMPTEL 1.8 MeV data. 

5. Conclusions 

An alternative imaging method for COMPTEL 7-ray 
data is presented, using a newly developed multircsolution 
reconstruction algorithm based on wavelets. The max- 
imum entropy and Richardson-Lucy algorithms, which 
have been used previously for COMPTEL image recon- 
struction, are very sensitive to statistical noise in the data, 
leading to image lumpiness and 'hot spots' in the recon- 
structed intensity maps. In particular, artificial 'hot spots' 
are indistinguishable from real point sources on basis of 
the sky maps alone, requiring an substantial additional 
analysis effort to assess their reality. We present the result- 
ing reconstructed image of the 1.809 MeV sky as an alter- 
native view, complementing previously presented images 
from the other methods, and pointing out their limita- 
tions. In particular, we caution to overinterpret structure 
in the ME and RL sky maps when modelling the galactic 
26 A1 emission for other studies (e.g. Lentz et al. 1998). 

Applying our new algorithm to COMPTEL data 
largely reduces or even removes artificial 'hot spots' and 
image lumpiness, depending on the selected significance 
requirement s. Simulations indicate that s = 3.0 seems to 
provide a reasonable choice for the reconstruction: while 
artifacts are mainly removed from the image, hints for 
weak (3 — 4cr) point sources are still present. Neverthe- 
less, it should be clear that the MREM sky map obtained 
in this work not necessarily provides a realistic view of 
the 1.809 MeV sky. The real 1.809 MeV intensity profile is 
probably much more confined to the galactic plane than 
the emission in the MREM sky map, but with an angu- 
lar resolution of 4° (FWHM), COMPTEL is not capa- 
ble of resolving this confinement. The 1.809 MeV emis- 
sion along the galactic plane may be much more struc- 
tured than shown in the MREM map, but the sensitivity 
of COMPTEL is not sufficient to map this structure. In 
this sense, MREM provides a more reliable image of the 
1.809 MeV 7-ray sky with respect to ME and RL since 
it does not show emission structures for which there is no 
strong evidence. Weak emission features which are close to 
the sensitivity limit of COMPTEL may however be sup- 
pressed in the MREM maps. Therefore ME and RL maps 
are used as complementary analysis tools. 
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